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An efficient metliod for calculating tlie electronic structure in large systems witli a fully converged 
BZ sampling is presented. The method is based on a k ■ p— like approximation developed in the 
framework of the density functional perturbation theory. The reliability and efficiency of the method 
are demostrated in test calculations on Ar and Si supercells. 



In periodic systems the correct description of electronic 
structure requires an exhaustive sampling of the Brillouin 
Zone (BZ). However for large systems, such as those that 
are employed in ab-initio molecular dynamics, this can 
be too expensive and very often the sampling is limited 
to the r point of the BZ. This is quantitatively accurate 
for very large systems and/or wide gap insulators. In 
all the other cases a correct sampling of the BZ is to 
be recommended. Another limitation of F calculations is 
the very coarse description that they give of the electronic 
density of states, even for rather large supercells. 

However, since full BZ sampling can be too costly, 
k ■ p— like approaches have been suggested by several au- 
thors in order to include k 7^ states at reduced com- 
putational costs ||l|(2|. In these works it has been sug- 
gested that the calculation of the k dispersion of eigen- 
values can be conducted in a restricted Hilbert space. 
Robertson and Payne |Q have proposed that this Hilbert 
space could be constructed from the occupied and a lim- 
ited number of unoccupied eigenfunctions of the F only 
KS Hamiltonian. Since by increasing the number of un- 
occupied states one converges to the exact result, this 
method is in principle exact. Unfortunately, the number 
of excited states that needs to be included in order to 
achieve a satisfactory convergence is large, therefore this 
approach has had limited impact. Kohanoff and Scan- 
dolo have proposed a procedure in which the basis 
vectors of the Hilbert space are variationally improved, 
along with the electronic wavefunction coefficients. Still, 
in the non-local PP formalism, this procedure would be- 
come rather expensive and actually it has been applied 
only to the local PP, which is quite restrictive. 

Here we adapt a recently developed version of density 
functional perturbation theory (DFPT) to cal- 

culate the finite k corrections. We show that this can 
give very accurate energies for systems of a few hundred 
atoms for a semiconductor like Si, and even smaller for 
an insulator like Ar, producing a significant improvement 
of the F only sampling at reduced computational cost. 
Furthermore we show that for large systems, the Hilbert 
space of the occupied eigenstates supplemented by the 



perturbative corrections provide an excellent description 
of the electronic density of states (DOS) at a much re- 
duced computational cost. 

We shall restrict ourselves to the case of a semiconduc- 
tor or insulator, where all states '!/'km(r) labeled by the 
state index m and by the BZ vector k are doubly occu- 
pied. We shall develop our method within the context 
of pseudopotential (PP) formalism, but it can easily be 
extended to other electronic structure schemes. In or- 
der to be definite we shall consider fully non-local PP 
of the Kleinman-Beylander type ||] composed of a local 
part, Vioc(r), and a non-local part, Y.i,l \Pi,l)^l{Pi,l\, 
where the sum runs over all nuclei / and the angular mo- 
mentum channels L. The Bloch theorem Q allows us to 
write V'km(r) = Mkm(r)e''' '' where Ukm(r) has the lattice 
periodicity. In terms of the Mkm the Kohn and Sham 
(KS) density functional ||^ can be rewritten as 

Eks = (1) 

BZ 

where the contribution from each k vector reads 
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The sum in m is over the occupied states. The electronic 
charge density is given by 

p(r) = 2E5]|^k™(r)P, (3) 

m BZ 

Ehxc is the sum of the Hartree and exchange and corre- 
lation contributions to the energy density functional. In 
the definition of the energy functional the orthonormality 



2 



condition is implied 



(4) 



For large unit cells, the BZ is very small and the func- 
tional in Eq. (^) can be expanded in k. The k dependence 
is in part explicit and in part implicit through the depen- 
dence of Ukm on k, which to linear order can be written 
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where the the 0*'' order wavefunctions, namely 

the results of a F only calculation, and the UamS the 
first order corrections. Since in the absence of a magnetic 
field, u^'m's and uim's can be made real the first order 
correction to the density vanishes identically. 
To the second order in k, £'(k) becomes: 
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Once the Urn S RrC calculated, the be deter- 

mined by minimizing the second order correction of the 
functional, E^'^\ relative to the Uam's. The minimiza- 
tion procedure is to be supplemented by the restricted 
orthonormality condition M 







y m,n and a 
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It can be seen that i?^^^ has the same structure as the 
standard form of the variational DFPT The only 
peculiarity is that the quadratic terms describing the ef- 
fect of the change in p due to the perturbation are ab- 
sent. Such a simplification is a consequence of the na- 
ture of the perturbation and occurs also, for instance, 
in the case of the calculation of NMR chemical shifts 
1^. In the DFPT the terms that are linear in {uaL} 
express the perturbative coupling m. In our case for a 



purely local PP the perturbation appears as a coupling 
to the current operator — iV, as is to be expected from 
a k ■ p scheme. Since it is well known that a non-local 
PP also affects the current operator, it is not surpris- 
ing that in the non-local case some extra perturbation 
terms appear. We note that as written in Eq.(^ the k 
expansion has only a formal character, since the posi- 
tion operator, in the non-local PP terms, is ill defined in 
a periodic system. If however we express the full non- 
local PP projectors in the reciprocal space and then ex- 
pand in k the non-local energy terms, we arrive at well- 
defined expressions that can be straightforwardly calcu- 
lated 1^ . In practice we first perform a F only calculation 
and then, using the variational DFPT module developed 
in our group ||5|, we perform three separate variational 
calculations to determine {u^m iU}. The calcula- 

tions of {wm^} and {uam} can take full advantage of the 
fact that these quantities are real. Once the {uim} are 
known, A^^ = 5^S/(i5fca5fc/3)|k=o can be evaluated and 
summing over the BZ we arrive at the following energy 
estimate 
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BZ Q/3 



Therefore, the computational effort required by the 
whole procedure is comparable to the optimization of 
four independent sets of wavefunctions jsj and it does 
not depend on the number of fc-points used in the BZ 
sampling. 

We tested the efficiency of the method on the crys- 
talline systems of Ar and Si. For these two systems we 
have considered supercells of increasing size and com- 
pared the F only calculation with and without the per- 
turbation theory correction to the exact result obtained 
performing a full BZ sampling in the primitive cell. 

The Ar fee unit cell, with 9.94 a.u. lattice constant, 
contains one atom and 8 valence electrons. We describe 



this system by means of a Goedecker PP |10 and the 
LDA functional. The converged total energy is achieved 
by the standard CPMD wavefunction optimization, 
over a (8 X 8 X 8) Monkhorst-Pack (MP) ||l|| mesh in 
the BZ of the primitive cell. Calculating the total energy 
with the F point only produces errors which decreases 
with the system size, as shown in Fig. |[ As expected in 
this case, since Ar bands are rather flat, the k-p approach 
converges very rapidly with the system size. 

We turn now to the study of the prototypical semicon- 
ductor. Si. The Si unit cell contains 2 atoms and 8 va- 
lence electrons. In this case we use the Troullier-Martins 
norm-conserving PP and the Ceperley-Alder local 
density approximation to the exchange-correlation term 
[Q. Convergence of total energy for this unit cell is 
achieved by an integration over the (15 x 15 x 15) MP 
mesh. As shown in Fig. ^, the DFPT correction produces 
an important improvement, such that the total energy er- 
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FIG. 1: Total energy vs the number of atoms in Ar super- 
cells. The energy differences are taken with respect to the 
exact result corresponding to the full BZ sampling. Circles 
represent differences due to T only calculations, squares show 
the effects of DFPT correction for a BZ sampling generated 
61/ (4 X 4 X 4) Monkhorst-Pack /uj/ mesh. In both calculations 
an energy cutoff of 80 Ry was used. 
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FIG. 2: Total energy vs the number of atoms in Si supercells. 
The tech nique details are the same as in Fig [J, but for the 
energy cutoff which is 20 Ry. 

ror is reduced by one order of magnitude for 128-atom 
size, and to less than 0.001% in the case of the 250-atoni 
supercell. Of course since the Si bands are much broader 
than those of Ar, convergence for Si is slower. 

Unfortunately the calculation of A^^ by DFPT pro- 
vides only the estimation of the total energy, which de- 
pends on the sum over the occupied state eigenvalues, 
and the contribution of the individual states cannot be 



separated. This is however necessary if one wishes to 
calculate the electronic density of states. We show here 
that a very good electronic density of states can be ob- 
tained if for any k one diagonalizes the k dependent KS 
Hamiltonian 

i7(k) = F(0) + ik^ - ik • V 

I,L 

- Y.\PilW^Pil\ (9) 

I,L 

in the Hilbert space spanned by the vectors {um''} and 
{"'^km}' which has dimension 2N, where 

"km = ik^u^xl + iky^^yl + ik^z^- (10) 

These vectors are orthogonal to the subspace {m™''} but 
they are not mutually orthogonal. We perform the or- 
thonormalization by means of the Lowdin method, which 
requires the overlap matrix 

to be calculated for each k vector. The matrices 
{u^am\u'^n) Can bc Calculated only once and then stored 
for later use. Similarly it can easily be shown that most 
of the matrix elements of i?(k) can be expressed in terms 
of k independent matrices, which can be calculated out- 
side the k loop. The only exception is the non-local 
PP term, which requires the evaluation of the projec- 
tors (P/L|e*"-|M^^) and {PiL\e^"'\ufj. However, fully 
in the spirit of our method, we can expand these projec- 
tors to second order in k, which reduces the evaluation 
to a linear combination of k independent quantities. 

We have tested our method on Ar and Si. Once again, 
for Ar the convergence is very fast and already for the 
32-atom supercell, the k-p approach gives results that are 
practically indistinguishable from the real ones. We illus- 
trate here the more interesting case of Si. In Fig. ||(a) we 
report the F point only DOS, as resulted for a 128-atom 
Si supercell. It is easily seen that this is a caricature 
of the real DOS. For the same system a much better de- 
scription can be obtained using a (4 x 4 x 4) MP mesh, as 
shown in Fig. |^(b). In the same picture we demonstrate 
that our method gives accurate results and that the use 
of the expansion in the evaluation of the non-local PP 
introduces only small errors. While the result of the full 
calculation and of the k • p method are comparable, the 
latter is computationally more efhcient by at least two or- 
ders of magnitude. This allows us to calculate easily the 
electronic DOS for the 128-atom supercell including 4096 
k points on an Intel Xeon PC. The result is compared in 
Fig.^ with the DOS obtained by a standard method on 
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FIG. 3: Electronic density of state calculated for the Si su- 
percell containing 128 atoms, a) V only DOS, h) Comparison 
between the standard CPMD (dashed) and the DFPT (solid) 
results, by using a (4 x 4 x 4) MP mesh. 



the primitive cell of two atoms and with 8000 equally dis- 
tributed k points. The only symmetry used is inversion. 
For both calculations the tetrahedron method was 
used in order to perform the k point integration in the 
BZ. The very good agreement between the two curves 
demonstrates the usefulness of our approach. 

In conclusion, we have shown that even for systems 
that are by present computational standards modest in 
size, the k • p method implemented here is accurate and 
computationally efficient. It is therefore now possible 
to perform efhcient ab-initio molecular dynamics calcu- 
lations with fully converged BZ sampling. 

We are very grateful to the Max Planck Institiit fur 
Festkorperforschung, (Stuttgart, Germany) for its sup- 
port during the initial phases of this research. 
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FIG. 4: Si electronic density of states: a) conventional fully 
converged calculation for a 2-atom unit cell with 8000 k- 
points, b) DFPT for a 128-atom supercell with 4096 k-points. 
In both cases the integration in the BZ by the tetrahedron 
method has been used. 
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